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[i] We present the results of a finite difference implementation of the kinetic 
Fokker-Planck model with an exact form of the nonlinear collisional operator. The model 
is lime dependent and three-dimensional; one spatial dimension and two in velocity space. 

The spatial dimension is aligned with the local magnetic field, and the velocity space is 
defined by the magnitude of die velocity and die cosine of pitch angle. An important 
new feature of model, the concept of integration along the particle trajectories, is discussed 
in detail. Integration along the trajectories combined with the operator time splitting 
technique results in a solution scheme which accurately accounts for both the fast 
convection of the particles along the magnetic field lines and relatively slow collisional 
process* We present several tests of the model's performance and also discuss simulation 
results of the evolution of the plasma distribution for realistic conditions in Earth’s 
plasmasphere under different scenarios. 
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I, Introduction 

[2] Magneto sphere -ionosphere (MI) coupling has interested 
scientists for decades, and in spite of experimental and theo- 
retical research efforts, ii remains one of the least understood 
dynamical processes in space plasma. The reason for this is 
that the numerous physical processes associated with MI 
coupling occur over multiple spatial and temporal scales. One 
typical example of Ml coupling is the production of up flowing 
ion events (or ionospheric outflows), such as auroral acceler- 
ation, ion energization in the cleft ion fountain, convective 
heating, polar wind, and plasmaspheric refilling. The classifi- 
cation of ionospheric outflows contributing to MI coupling 
has frequently been divided into two broad physical catego- 
ries, The first of these is the polar wind which exists on high- 
latitude field lines connecting to the interplanetary medium 
and the geomagnetic tail* Several tutorials by Schunk [1986, 
1988a, 1988b] and reviews by Ganguli [1996] and Yau et al 
[2007] provide a complete picture of the historical develop- 
ment of polar wind studies. The second region where iono- 
spheric outflows are important is the plasmasphere, where 
closed field Sines allow' this region of the inner magnetosphere 

to become saturated with thermal ions (e.g., see the review by 
Singh and Horwitz [ 1 992]). 
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[3] The present generation of ion up flow models is based 
either on a truncated series of higher-order velocity moments 
[e.g., Schunk and Watkins, 1979; Barakat and Schunk, 
1982a, 1982b; Mitchell and Paimadesso, 1983; Khazanov 
et ah, 1984; Gombosi et «/ , 1985; Singh and Schunk, 1985; 
Combos i and Killeen , 1987; Ganguli and Mitchell , 1987; 
Ganguli et al, 1993; Rasmussen and Schunk, 1988; Demurs 
and Schunk, 1 989, 1 994; Moffett et al, 1 989; Singh and Torr, 
1990; Kdtvsmezey et al , 1992] or on kinetic methods 
including simplified hybrid PIC simulations [e*g., Barakat 
and Lemaire , 1990; Wilson et ai , 1990; Demurs and 
Schunk, 1989, 1992, Wilson , 1992; Barakat et al ., 1993, 
1995; Miller et al . , 1993, 1995; Barakat and Schunk, 2001; 
Horwitz and Zeng, 2009; Barghouthi et al . , 201 1] and direct 
solution of die kinetic equations [e.g*, Lemaire and Scherer, 
1972; Lie-Svendsen and Rees , 1996; Khazanov et a!., 1997]. 

[4] These techniques have been further applied to study 
the global nature of ion up flow [e.g., Schunk and Sojka, 
1989; Gardner and Schunk* 2004; Barakat and Schunk, 
2006; Glocer et at., 2009]. Both of these methods have 
powerful strengths and considerable weaknesses dial have 
been reviewed and discussed by Echim et a L [2011]* 

[5] The higher-order fluid models have serious limitations 
when they are apphed to regions where collisions are infre- 
quent or negligible. These limitations are a result of a fun- 
damental approximation used by all generalized transport 
models; the models are based on a perturbation approach 
which assumes distribution functions are close to Maxwellian 
or bi-Maxwdliam This approximation makes fluid models 
computationally efficient, but limits their applicability at 
higher altitudes* The bi-Max wellian based perturbation 
approaches are more appropriate for gyration dominated 
plasma, but still require collisional dominance* It is well 
knowit dull collisions drive the distribution function toward 
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equilibrium: this process is the physical reason behind all 
velocity moment based approximations. A natural consequence 
of this is that as collisions become less and less frequent, the 
velocity distribution can develop highly nonequilibrium fea- 
tures (such as conics or double humps) that cannot, be accounted 
for with perturbation methods. To put it plainly, generalized 
transport methods lose validity in collisionless regimes and 
must be replaced by a kinetic treatment. 

[6] True kinetic methods provide a full solution for the 
multispecies phase- space distribution function with respect 
to seven independent variables: time, three-dimensional 
location, and the particle velocity vector. Results from a 
kinetic solution can be directly compared with in situ 
observations of the ion distributions by spectrometers on 
rockets and satellites with no need to take moments of the 
measurements. The full informational content of the data 
may be exploited in this type of comparison. By looking for 
characteristic “signatures” in the distribution function, one 
would he able to identify which physical mechanisms are 
responsible for certain ion outflow events, 

[?] The disadvantage of a kinetic treatment is that it is 
computationally demanding, and that development of ade- 
quate solution methods is not straightforward. Presently no 
existing plasma outflow model is based on such an approach 
and it is not likely that anyone will be able to develop one in 
the near future. Also, most of the previous kinetic iono- 
spheric outflow studies have used the statistically based 
Monte Carlo technique [e.g., Bamkat and Lemaire. 1990; 
Wilson , 1992; Miller et al., 1993; Barghouthi et al , 1993], 
This method is encumbered with random uncertainties from 
the particle simulation, as described by Miller and Combi 
[1994] and Bamkat et al. [1998]. Therefore, great care 
must be taken when applying this technique, 

[s] Another approach, and the one applied in this study, is 
to use the Fokker-Planck collisions! operator, similar to the 
studies by Khazanov et al. [1994] and Lie-Svendsen and 
Rees [1996], The main development of this paper and the 
difference from all previous studies, however, is the use of 
the exact fomi of the Fokker-Planck collisional operator 
without any assumptions with respect to the distribution 
function of the background particles. In this study the non- 
linear form of the Fokker-Planck collisional operator has 
been solved self-consistent ly for the first time for space 
plasma applications. This method will not have to deal with 
the statistical uncertainties of a particle simulation because 
the solution procedure of such a problem is not statistical in 
nature [see Khazanov et al ,, 1994; Khazanov and Liemohn, 
1995]. For the sake of simplicity, this paper will only 
focus on the electron distribution function formation, leaving 
the discussion of electrodynamic and Coulomb electron-ion 
coupling to our forthcoming paper. It should be noted 
though that the numerical approach developed here will be 
exactly the same for the ion population of space plasma and 
is therefore presented here in rather general form. 

2, The Fokker-Planck Kinetic Equation 

[9] As mentioned above, we will use the exact form of the 
Fokker-Planck collisional operator to develop a complete 
kinetic description of ionospheric plasma outflows from the 
collision-dominated region to the collision less magnelo spheric 
plasma. Such an approach provides a continuous calculation of 


the self-consistent coupling processes between the different 
components of the ionospheric and magnetospheric plasmas 
along geomagnetic field lines. The Fokker-Planck kinetic 
equation, written in terms of the two Rosenbluth potentials 
HJy) and G ft (v), in the presence of gravitational g, electric E 
and magnetic B fields can be presented in the following gen- 
eral form [Rosenbluth et al . , 1957; Shkarofsky et aL 1966]: 


g,vv r /t[, + i(EAxs)].v,/-£: 
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I lere/=/(A r, v) is the distribution function, where r and v are 
the position and velocity vectors, respectively; index a 
denotes the background species with which the particle with 
charge e and mass m collides; and n a is the density of the 
species a. The Rosenbluth potentials H a (v) and G r> (v) are 
integrals over the distribution of the background particles 


^ f /.( v ')l v ~ VV 3 ' 


( 2 ) 


[ 10 ] Previous use of the Fokker-Planck collisional operator 
for ionospheric plasma outflows was restricted to its simpli- 
fied, linearized, form by assuming that the thermal background 
electrons and ions are static and have Maxwellian distribution 
functions. This greatly simplifies the calculation of the 
Rosenbluth potentials and equation (2) can be presented in an 
analytical form. Such a linearized Coulomb collisional oper- 
ator has been used in earlier calculations of superlhemial 
electron transport [Khazanov et al y 1 979; Yasseen et at., 1989; 
Khazanov et ai f 1994] and in the polar wind model by Lie - 
Svmdsen and Rees [1996]. It should be stressed, however, 
that in some eases the departure of the plasma distribution 
function from a Maxwellian or bi-Maxwellian cart be very 
large, causing the linearized Fokker-Planck collisional opera- 
tor to lose validity. For example, [Bara hat et al. , 1995], using 
their Monte Carlo model, found that between the low-altitude 
collision-dominated and high-altitude collision less regions, 
the If velocity distribution becomes double-humped in 
energy. The formation of this double hump is a natural con- 
sequence of the interplay between the electrostatic ion accel- 
eration and the velocity- dependent Coulomb (IT -CT) 
collisions. It may also occur in other regions of space plasma 
and should be rigorously analyzed with the model being pro- 
posed in this study. This will be the case not only for the ions 
but also for the electrons. Such a unified treatment for the 
electron distribution function is especially needed in the 
presence of superthermal electrons. Artificial separation of 
the electron distribution function into thermal and super- 
thermal parts Leads, in some cases, to unrealistically high 
values of electron temperatures and panicle fluxes [Tam etal 
1995], 

[n] As we pointed out in section 1, we use the exact form 
of the Coulomb collisional operator (!) without any 
assumptions with respea to the distribution function of the 
background particles. In this case, the collisional operator 
becomes nonlinear and depends on the plasma distribution 
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Figure L Geometry of Lhe problem. Spatial coordinate s 
along the magnetic Held line. Principal variables in velocity 
space are the velocity of the particle v and cosine of pitch 
angle £ = cos 9. 


function itself because collisions between similar particles 
may no longer be neglected. Also, we expect that in some 
cases the departure of the plasma distribution function from 
Maxwellian can be very large and the linearized Fokker- 
Planck colhsional operator can lose validity, even for the 
interaction between different kinds of particles. To deal with 
this problem, we will transform equation ( l ) into a computa- 
tionally manageable form below, as established and described 
by Khazanov et at [1979, 1994, 1996; Khazanov and 
Liemohn, 1995] for the case of superthermal electron trans- 
port. Nonlinearity in equation (I) will be handled using an 
iterative scheme similar to the case w hen we used the isotropic 
part of the Fokker- Planck collistonal operator to calculate the 
electron distribution function in the collision-dominated 
region [Khazanov et at., 1978. 1979]. 

[ ] 2 .] Another type of nonlinearity occurs in our model 
through the development of the self-consistent electric field 
that is part of the calculation of the hydrodynamic model and 
the Fokker-PIanck kinetic equation (i). The reason for the 
formation of a self-consistent potential in a collisionless 
plasma is this: high mobility electrons lend to overtake ions. 
As a result, the electric neutrality of the plasma is violated 
and an electric Held appears which constrains the electrons, 
forcing them, on average, to travel together with the ions. 
This field also significantly affects the motion of the ions by 
accelerating them. The electric field acts as a catalyst by 
transferring Lhe pressure of the electron gas to the plasma 
ions; this pressure is proportional to the electron temperature, 
7 1 ,.. Therefore, when T e - T u the effects of the self-consistent 
electric field and the effects of tlie ions 1 thermal motion are 
generally of the same order of magnitude. Photoelectrons, 
which form due to ionization of the atmosphere by solar 
radiation, can alter the self-consistent potential in the space 
plasma. The presence of the enhanced high-vdocity tail in 
lhe electron distribution will increase the number of fast ions. 
Due to the enhanced ion acceleration in an expanding 
plasma, the initial superthermal electron distribution function 
could be changed. As we pointed out in the Introduction, the 
electrodynamic electron- ion coupling will be ignored in our 
current study for sake of simplicity. The calculation of the 
sell -consistent electric field wall be included in the iteration 
loop of the model in future studies, similar to the descriptions 


of Khazanov ei ai [ 1 997], based on the quasi -neutrality and 
current less conditions. 


3. Description of the Model 

[□] In the Earth’s strong magnetic field the distribution 
function is highly symmetric in the plane transverse to the 
magnetic field direction. It is therefore convenient to choose the 
spatial coordinate V along the local magnetic field line and in 
the velocity space to choose the spherical system of coordinates 
with polar axis along the local magnetic field line (Figure l). 

[\a\ Under assumption of azimuthal symmetry of the dis- 
tribution function the equation (t) is Iran s formed to the form 


df , ¥ 

— -E-fv— — 
dt 5 ds 


\ I dB df 

2 V 5"&3f 


= Hf) 


(3) 


Here we omitted gravitational and electric fields. 'Hie vari- 
able £, is Lhe cosine of the angle between particle velocity 
vector and magnetic field direction, or in other words it is 
cosine of the pitch angle, which in our geometry corresponds 
to the polar angle. Only in homogeneity of the magnetic field 
is included in the last term of the left hand side. This term is 
responsible for conservation of Lhe adiabatic moment of the 
particle in the absence of collisions 

= const (4) 


and represents the most important dynamic effect of trapping 
particles in the Eartlris magnetic field. 

[is] The Coulomb colhsional operator will be in its exact 
form, similar to [Khahihrakhmanov and Khazanov, 2000], 
which in spherical system of reference in velocity space 
takes the following form: 
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Here, colhsional strength is defined by 
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and Fokker-PIanck coefficients can be expressed as 
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in terms of the Rosenbluth potentials G and H. The integral 
definitions (equation (2)} for these potentials are equivalent 
to couple of Poisson equations in the velocity space 

A V H = -/, A ,C = // (12) 

As was originally suggested by Rosenbluth et ai [1957], it is 
convenient to use a truncated series of Legendre poly- 
normals L*{0 for the expansion in the angle variable £ 

/(v,tf.= 5>(v)4(0 (13) 

A-=0 

This type of expansion has been widely used in fusion 
research [e.g., Killeen et al , 1976], where magnetic particle 
trapping and collisional transport is very sensitive to fine 
details of the plasma particle distribution function* Khazanov 
[1979] also used this approach in order to describe ionosphere- 
plasmasphere transport of superthermal electrons* Pierrard 
and Lemalre [1998] as well as Pierrard et al [2001] used 
similar expansions to study the polar wind and, more recently, 
to obtain a self-consistent description of the electrons in the 
solar w ind. 

[i6] For the expansion coefficients H k {v) and G*(v) 
then we have a set of second-order ordinary differential 
equations: 


SSE 

V 

1 

II 

£ 

1 

1 

rt >i<§ 

(14) 

a a 

— V 2 — G k - k(k - 1)G* = \rH k 
&v 8v ’ 
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with appropriate boundary conditions. The solution again 
can be expressed in terms of corresponding Green’s function 
[Rosenbluth et. ai, 1957], for numerical purposes; however, 
it is more convenient to solve corresponding boundary value 
problems using standard finite difference approximation* 

[i 7] It is convenient at this point to transform to dimen- 
sionless variables. The spatial variable is normalized to the 
length of the given magnetic field line S so that dimension- 
less length along the field line varies from -I to l when 
particle moves from the boundary ia Northern Hemisphere to 
the conjugate point in the Southern Hemisphere 

s — ► — 1 + 2s l S 

We fix the maximum value of the velocity of the particle on 
the grid V With this choice the characteristic timescale of the 
problem becomes S / V, which is the time required for fastest 
particles under consideration to travel the distance S. In 
dimensionless variables equation (3) still has the same form 
with dimensionless collisional strength defined as 

47rc 4 5h!nA ,, 

r= (](>) 

Here we also assumed that distribution function is normal- 
ized to unity. Therefore, the problem is characterized by only 
one dimensionless parameter, the relative strength of the 
collisional operator. This is very important step. As a result, 
every numerical solution presented in the paper represents a 
whole family of solutions. Changes in physical plasma 


density n and scaling velocity V such that n / f 4 leaves the 
parameter ! ' unchanged do not affect the solution of the 
problem. Of course the interpretation of the result will change 
as the "kinematic" timescale is still defined by r = S / V and 
the physically most meaningful interpretation of the dimen- 
sionless parameter T is the ratio of that “kinematic” Litnescale 
to collisional timescale. When F < 1 is small w r e expect most 
of the particles collide very rarely with distribution function 
deviating strongly from Maxwellian. In opposite regime T 
l the distribution function is expected to stay close to Max- 
wellian form and the details of the deviation of the distribu- 
tion function are folly described by the value of the parameter 
and physical nature of the operator in the right hand site. The 
results are also invariant under change of foe particles mass 
and charge which conserve the factor IT 
[is] Equation (3) is subject to boundary conditions. In the 
spatial variable we specify the distribution function at s = - I 
and s = L In the velocity space it is required that the flux of 
the particles across the boundaries vanish. For the angular 
variable corresponding duxes vanish automatically as the 
F okker- P lanck c oeffic i e nts EL, are equal zero there. 

In velocity variable some of the Fokker-Planck coefficients 
are singular at v - 0 S for instance D ^ and ihe boundary 
condition 


must be imposed. In the past, Khazanov [ 1 979] and Pierrard 
and Lemalre [ 1998] used this “regularity condition" to avoid 
singularities which appear at v = 0, 


4. Numerical Implementation 

[ 1 9 ] N umerica 1 implement at i on o f the eq uat ion (3 ) fo 1 lo ws 
the general scheme of time splitting for multidimensional 
problems [Marchuk, 1975]. 

4J. Integration Along the Trajectories 

[ 20 ] In the absence of collisions the left hand side of 
equation (3) conserves the adiabatic moment of the particle 
given by (4). Therefore, for a given magnetic field profile B 
(s) the particle moves in the plane (s\ £) along the trajectories 
defined by the integral of motion (4), These trajectories are 
independent of particle velocity v and are shown in Figure 2 
for dipole magnetic field as the lines connecting die grid 
points. If we choose particular grid point { 5 ; A) in this plane as 
shown in Figure 2, the points (s T I ; k + 1 ) and (s — I ; k - 1) 
belong to the same trajectory. We can expand the value of the 
distribution function /,. + 1; A + j at the point (s T I ; k + 1) in 
Taylor series 



fyj-H Jf-l) 

$ik 


— £*_i) + ... 


(17) 


From expression (17) the centered spatial derivative in the 
(s, () plane al the grid point (.r; k) is approximated as 


r f , /i+i.*+i ~/.-u i 


'£\ kti - &- 1 
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Trajectories of the particle 
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Figure 2, Trajectories of the particles in (s, 0 plane as defined by the conservation of the magnetic 
moment of the particle. 


here we also look into account that £ = £($) along the trajec- 
tory and used new variable z 



Using approximation ( 1 8) in the left hand side of (3) 


(|) 



3A 

-&-i 
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_ Zf+] 
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3 fA,* z ‘+i 



\-£ I dB df\ 

2 B ds dU, k 


(19) 


one can see that the last two terms cancel each other when 


(i-f i'«') 

I I ”“ 41-1 \ 2 ^ ^ / sjt 


(20) 


But this is just a finite difference approximation of the 
equation for trajectories (4). In other words, when grid points 
(jt i , k 1 ) and (s + 1 , k + 1 ) belong to the same trajectory 
as point (s, k ), the left-hand side of equation (3) has simple 
finite difference approximation 



and conservation of the adiabatic moment of the particle is 
accounted by nonuni form grid, where all the grid points 
belong to trajectories of the particle. As we show section 4.2, 
this procedure is able to completely eliminate numerical dif- 
fusion in the absence of collisions. We note that Lie-Svendsen 
and Rees [1996] also obtained a solution to the Fokker- 
Planck equation for the case of polar wind minor ion outflow* 
but they used a standard finite difference rather than our 
approach which chooses grid points belonging to the 


trajectories of the particles in absence of collisions. More- 
over, their approach was for a steady state solution, while our 
approach is for a time-dependent solution. 

42. Time Splitting 

[ 21 1 The time splitting technique allows one to use inde- 
pendent approximations of the differential operators acting 
in each of the independent variables. Equation (3) is repre- 
sented as 


¥ 

3i 


M 

=EV 


Then M equations 


¥ 

dt 


= LJ f 


(22) 


(23) 


are solved in succession using as an initial value result of the 
previous fractional time step /J_ t . This technique is very 
flexible in allowing addition of new terms. The second -order 
approximation in time can be achieved by proper choice of 
the operator splitting in representation (22). 

43. Testing the Code 

[ 22 ] In this section we will validate the quality of the 
numerical approach. In particular, we will carry out four 
carefully crafted tests that demonstrate the model's fidelity 
to the well known properties of the Coulomb collisional 
operator, the model’s low numerical diffusion, and a com- 
parison with previous models. The following list provides a 
short description of each test with the complete results and 
details left to sections 4.3.1-43.4 

[23] 1 . One of the most important properties of the Cou- 
lomb operator is conservation of density and energy of the 
distribution during the process of relaxation toward a Max- 
wellian Ibrm. We test the Coulomb collision operator with a 
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t=1.0 


t=4.0 



Figure 3. Relaxation of the initial streaming distribution function (24) toward Maxwellian. 


fairly arbitrary initial distribution and follow its relaxation to 
a Maxwellian distribution (section 4,3.1), 

[ 24 ] 2, An important property of the transport operator in 
the left hand side of equation (3) is the complete isolation of 
i he trapping region from the loss cone. In the absence of the 
collisions, particles are trapped in the equatorial region or are 
freely moving along a magnetic field line (tom one end to the 
other; there is no diffusion into or out of the trapped region. 
We test this property by turning off the collision terms and 
considering an upflowing population in the ionosphere that 
completely tills the loss cone region (section 4,3.2). 

[is] 3. Inside the trapped region any initial nonuni formity 
of the distribution function along the trajectory of the parti* 
cie becomes quickly ‘‘bounce 1 ' averaged and uniform. We 
illustrate that important property by considering convective 
phase mixing of the distribution in the absence of collisions 
(section 4,3.3). 

[ 26 ] 4. Sharp boundaries can appear between trapped and 
free motion regions as result of collisions. This effect was 
the subject of numerous studies in the past. We compare our 
model results against these previous studies as an additional 
test (section 4.3*4). 

4.3.1. Relaxation to Maxwellian Distribution 

[ 27 ] In this test the convection term was turned off. It is 
assumed that distribution function is independent of the 


spatial variable but has an arbitrary dependence on velocity 
and pitch angle. Figure 3 consists of several snapshots of the 
time evolution of an initial distribution de lined by 

^(£,v) = IQv 2 cxp(-yVO.Q8k 2 (24) 

The space-pilch angle grid corresponds to L - 2 with total of 
75 grid points at the equator. Tire number of velocity grid 
points has been taken to be 64. Although the problem has 
been solved in spherical system of coordinates in velocity 
space, the results have been transformed lo v ± = vcosfl, vy = 
vsin# space for presentation. In this test, special attention has 
been drawn to the conservation properties of the numerical 
model. It was found that density and energy of the distri- 
bution is conserved within a fraction of a percent on collt- 
sional Lime scales (Figure 3), 

[zs] The final distribution is shown in Figure 4. Here the 
difference between calculated distribution function and 
Maxwellian 0.9 expf-v 2 / 0.08) has been plotted for every 
value of the pitch angle. As we can see, the distribution 
function has relaxed to a Maxwellian with relative error of a 
fraction of a percent for all pitch angles. 

4,3.2. Loss Cone Passing With No Collisions 

[ 2 ^] The second new feature of the model is integration 
along the trajectories. This technique helps to virtually 
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Difference with maxwellian 



Velocity, v 

Figure 4. Relaxation of the initial distribution function: 
approximation error. 


eliminate the numerical diffusion across the trajectories in 
the absence of interparticle collisions. This important prop- 
city is illustrated in Figure 5, Here the collisions have been 
turned off, and we start the calculations with an empty mag- 
netic field line. At one of the ionospheric boundary the dis- 
tribution function is specified such that the loss cone is 
populated but not the trapped region. As is clear from the 
results, the trapped region remains completely empty every- 
where along the Held line with a step-like increase at (he loss 
cone boundary, 

4.3*3* Convective Phase Mixing of the Distribution 
in the Absence of Collisions 

[30] As the second test of die concept of integration along 
particle trajectories the evolution of the initial distribution 
function with the maxima inside the trapped region has been 
modeled. Figure 6 (top left) shows an initial distribution 
which is localized at equator inside the trapped region and is 


moving toward the right boundary. The distribution is con- 
verted toward the right boundary as can be seen in Figure 6 
(top right). Later on, particles are reflected and move back 
toward the left boundary, as Figure 6 (bottom left). The 
circulation along the trajectories inside the trapped region is 
differential, Particles with larger velocity move faster along the 
trajectories and the initially localized maxima, even for the 
particles with the same velocity v as in Figure 4, is being dis- 
persed by the differential circulation. At later times (Figure 6, 
bottom right), the distribution attains a characteristic "‘hat-like'’ 
shape, where distribution is completely smoothed along the 
particle trajectory'. As before, in the absence of collisions there 
is virtually no diffusion of the particles across the boundary of 
the trapped region, 

[31] The effect of rapid phase mixing of the particles 
inside the trapped region is the basis for the so-called 
* "bounce-averaged 1 * description of the plasma, where the 
assumption that the distribution function is constant along 
ihe trajectory of the particle in s space considerably sim- 
plifies the description. 

4.3*4. Comparison With Earlier Model 

[ 32 ] Khazanov et al [1993] used the linearized coulomb 
collisions] operator in order to study the evolution of 
superthermal electron component. In their model the mir- 
roring term in equation (3) has been eliminated by direct 
change oT angle variable £ lo £0 according to the adiabatic 
invariant (4), 

[ 33 ] In order to verify the concept of integration along the 
trajectories the test calculation was performed on exactly 
same, linearized model of Khazanov et al. [1993, equation (1 )] 
but using algorithm presented in this paper. The result of the 
comparison is shown in Figure 7 as an equatorial distribution 
function for high-energy 50 eV electrons against the pitch 
angle at the L shell L - 3, As shown in Figure 7, integration 
along the particle trajectories yields a result very close lo dis- 
tribution obtained earlier by a different technique. 

5* Results and Discussions 

[34] The above tests of the numerical model has demon- 
strated that it is capable of accounting for particle convection 


a) 


T = ZQ r V=0,5 


b) T = 690, V=0,5 



Figure 5. Refilling of the empty held line tube in the absence of collisions. There is no numerical di ffu- 
sion Into the trapped region. 
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Figure 6. Phase mixing of the initial distribution function with maximum at the equator. 


without introducing spurious numerical diffusion between 
the loss cone and trapped region. In this section we present 
results of several simulations of the dynamics of the particle 
distribution function in the plasmasphere. We assume a 
magnetic dipole configuration and L = 2. In general, how- 
ever, this model can be adjusted to use an arbitrary magnetic 
held configuration* 

5.1. One-Sided Refilling of the Initially Empty 
Plasmasphere 

[ 35 ] As the plasma particles are treated in our model in a 
self-consistent and unified fashion it is now possible to 
compute the refilling of completely empty plasmasphere 
with a given source at the ionospheric boundary. The models 
with linearized version of collisional operator always assume 
the presence of a background plasma and therefore are nol 
capable of modeling a completely empty plasmasphere. 
Moreover, such models are incapable of handling the refill- 
ing of the core low-energy plasma* 

[36] For L = 2 and an initially empty magnetic 11 ux tube, 
the distribution function at the ionospheric boundary is taken 
in the form 

0.1 exp(-v*/0.05) (25) 


There is no specification of the particles in our model, which 
is described by only one dimensionless parameter (16). 
Thus, as we mentioned above, the developed approach is 
applicable for electron and ions. In Figures 8 and 9 the time 
evolution of one-sided refilling for parameter value T = 1 


Comparison 



Pitch angle, 6 

Figure 7* Comparison of the stationary distribution lunc- 
lion for 50 eV electrons at the equator /, = 3 with that com- 
puted by Kkazanov et a i [1993]. 
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T = 2.0 


T = 4.4 



T s= 20.4 


T = 40.0 



Figure 8. One-sided refilling of the empty plasmaspheric tube with L 2; distribution function at v 0.5 
as function of the distance s along the field line and cosine of the pitch angle 


is presented. Figure 8 (top left) shows the earlier stage of the 
refilling process. The loss cone is quickly filled by the 
source at the left boundary, while the trapped region remains 
essentially empty. At later times the small amount of slow 
particles, which has entered trapped region due to finite 
coilisional strength, appears as a component moving back* 
ward, from the right to the left (negative £). Particles accu* 
mulaie in the trapped region, with the peak intensity at the 
boundary of the trapped region, 

[ 37 ] Finally, a quasi -stationary state is approached, where 
the detailed balance between convection and coilisional 
diffusion in velocity space gives rise to the distribution 
function, which is far from a Maxwellian distribution. This 
is oven more clear from the high energy part of the distri- 
bution function shown in Figure 10 

5*2, Relaxation of an Initially Maxwellian Distribution 
in the Plasma sp here 

[va] In order to understand the importance of convection 
and particle trapping, the model with an initial isotropic 
Maxwellian distribution has been calculated. The distribu- 
tion at the boundary and initial distribution throughout the 
magnetic tube line is Maxwellian, The density of the initial 
distribution is nonuniform along the field line and is pro- 
portional to magnetic Held strength B(s)< Evolution of ihe 
distribution is shown in the same format as before. Figure 1 1 


shows the distribution of particles with v = 0.67 for four 
different times. The distribution near the right boundary of 
the magnetic tube line, s = —0.96 is shown in Figure 12, The 
redistribution of the particles during the relaxation is clear 
with a transient appearance of anisotropy and overall change 
in shape. Although the distribution for long time evolution 
appears to be isotropic at tow and intermediate energies, the 
strong anisotropy builds at high energies. This is shown in 
Figure 13 where the high-energy part of the same distribu- 
tion function as is shown in Figure 1 1 (bottom right). 

5*3* Particle Precipitation Into the Loss Cone With an 
Initially Maxwellian Distribution Function of Particles 

[ 3 y] The process of particle precipitation from the plas- 
maspheric tube is modeled again with an initially isotropic 
Maxwellian distribution everywhere in the tube. The evo- 
lution of the computed distribution function at velocity v = 
0-67 is shown in Figure 14 r The particles quickly precipitate 
from the loss cone. It takes much longer for particles in the 
trapped region as a result of collisions to scatter to loss cone 
and precipitate. The velocity distribution function dose to 
the right boundary at s = -0,96 is shown in Figure 15. It is 
dear that a strongly anisotropic distribution is formed due to 
the precipitation near the ionospheric boundary. The high- 
energy tail of the distribution at T = 800 is shown separately 
in Figure 16 and demon si rates that the anisotropy of the 
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Figure 9, One-sided refilling of the empty plasmaspheric lube with L ~ 2 distribution function at s = 
-0.96 as function of the velocity v and cosine of the pitch angle £ 


Equatorial DtstribLiEon, T - 40,0 



Figure 10, One-sided refilling of the empty plasmaspheric tube with L — 2: high-energy part of the dis- 
tribution function at s - -0.96 as function of the velocity v and cosine of the pitch angle £. 
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Figure 1 1 , Relaxation of the initial Maxwellian distribution function to stationary distribution at v = 0,67 
as a function of distance s and the cosine of the pilch angle £ = cos 0. 
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Figure 12, Relaxation oF the initial Maxwellian distribution Function to stationary distribution at s 
—0.96 as junction oF velocity v and the cosine of the pitch angle £ = cos 0. 


Distribution at s ^ -0.96, T = 400 



Figure 13* I ligh-energy pari of the stationary distribution at s = 0.96 as function of velocity v and cosine 
oF pitch angle £ = cos 0. 
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Figure 14* Precipitation from the plasmaspheric tube with L = 2: distribution function at v = 0,67 as a 
function of distance s and cosine of pitch angle £ = cos 0, 
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Figure 15* Precipitation from the piasmaspheric tube with L 2: distribution function at s = -0.96 as a 
function of velocity v and cosine of pitch angle £ — cos 0. 
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distribution function is very strong due to the very slow 
collision rate of high-energy particles. 

6. Con elusions 

[ 40 ] The kinetic Fokker-Planck model with an exact non- 
linear Coulomb collisional operator has been implemented 
numerically and tested. The concept of integration along the 
particle trajectories and time splitting technique has been 
developed and used. This allows us to obtain an accurate 
description of collision a! relaxation of the particles in the 
presence of fast convection along the field lines in the loss 
cone, rapid circulation inside the trapped region, and relatively 
slower collisional diffusion of velocity space. The numerical 
model has been used to describe the plasma dynamics and 
transport in the plasmasphere between the conjugate regions of 
the ionosphere. The numerical modeling results have been 
presented in the form of time-dependent distribution functions 
of the charged particles along the geomagnetic dipole mag- 
netic Held line as a function of velocity and pitch angle. 
Results fur high -energy tails of the distribution are compatible 
with the previous models, based on a simplified description of 
the Coulomb collisions. 

|4i] The use of exact form of Coulomb collisional opera- 
tor allows us to overcome most serious limitations of the 
linearization procedure used in previous research. As a 
result, a complete picture of particle distribution function is 
available. This represents a major step toward the accurate 
self-consistent global modeling of the ionospheric outflows 
with inclusion of self-consistent electromagnetic effects and 
multiple plasma components. 
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